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PREFACE 


It is well known that a calculus-of-variations approach to 
solving the Bolza form of trajectory optimization problems usually yields 
a nonlinear two-point boundary value problem in terms of the state and 
Lagrange multiplier variables. Closed-form solutions for problems of 
this type are difficult to obtain except for a few simple problems. As 
a result, recent work in trajectory optimization has focused on numerical 
procedures for obtaining solutions using high-speed digital computers. 
Particular interest has centered on a group known as the second-order 
methods . 

One such method is the Successive Sweep Method (SSM). It uses 
the generalized Riccati transformation technique to bypass a direct nu- 
merical integration of the perturbation equations. The reasons such an 
approach has much potential appeal are presented in this study; however, 
because the SSM iterates on the control values over the interval of in- 
terest, considerable computer storage is necessary even for problems of 
small dimension. This storage is required to compute corrections to the 
assumed control. Furthermore, the Eulerian control is not obtained upon 
convergence . 

This research develops a new second-order numerical optimization 
method, the Modified Sweep Method (MSM) . It requires very little computer 
storage and provides the Eulerian control. In addition, the properties 
and information contained in the Riccati transformation variables are 
preserved. 

Furthermore, this research also presents a new scheme for defi- 
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ning classes of numerical optimization methods. The Successive Sweep 
Method and the Modified Sweep Method are then discussed in terms of 
differences arising because each falls into a different class. The Modi- 
fied Sweep Method is subsequently compared numerically to the Method of 
Perturbation Functions (MPF) , both of which belong to the same class. 

The author extends thanks to Mr. Walt Williamson of The Univer- 
sity of Texas at Austin for many helpful discussions concerning the Apollo 
reentry problem, to Mr. I. J. Kim, of Lockheed Electronics Company, Hous- 
ton, for programming the plots, and both to Mr. Kim and Mr. Mike Frederick 
of The University of Texas at Austin for helping with the data. He also 
expresses gratitude to Dr. J. M. Lewallen of NASA/MSC and Professors A.M. 
Bedford and E.J. Prouse of The University of Texas at Austin for serving 
on the dissertation committee and helping with the manuscript. Special 
thanks is due to Mr. E. L. Davis of NASA/MSC for making this research 
possible, and for his personal interest and friendship which have pro- 
vided constant inspiration. 

The author expresses his indebtedness to Professor B. D. Tapley 
of The University of Texas at Austin, who suggested this research and 
provided many stimulating discussions while serving as committee super- 
visor. 
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A MODIFIED SWEEP METHOD FOR 


CONTROL OPTIMIZATION 


Daniel Colunga, Ph.D. 

The University of Texas at Austin, 1970 

Supervising Professor: B. D. Tapley 

A Modified Successive Sweep Method was devised which yields 
Eulerian solutions to two-point boundary value problems of control opti- 
mization. This was accomplished by requiring control satisfaction of 
local optimality over the entire time interval of interest while simul- 
taneously relaxing terminal transversal! ty requirements on the Lagrange 
multipliers. 

The new method was tested successfully on several classes of 
problems including optimizing the roll program for an Apollo-type three- 
dimensional reentry trajectory so as to minimize a time integral of 
spacecraft heating and acceleration. 

This new method was shown to require significantly less computer 
storage than the original Successive Sweep Method while requiring numeri- 
cal integration of fewer variables. In addition, the method was shown to 
possess rapid terminal convergence and a conjugate-point test capability. 
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CHAPTER 1 


INTRODUCTION 

Optimal control of spacecraft trajectories requires obtaining 
an optimal (maximum or minimum) value for an appropriate scalar quantity. 
This scalar quantity measures spacecraft performance and is called the 
performance index. In addition, terminal conditions such as might be 
specified for intercept or rendezvous problems must often be satisfied 
simultaneously. 

After the problem has been formulated mathematically, several 
conceptual approaches are available to obtain the conditions required to 
solve the optimization problem. Among the more usual approaches are the 
calculus-of-variations , dynamic programming, and Pontryagin's Principle. 

The calculus-of-variations is considered here because specific 
optimal control problems can be considered as particular cases of the 
more generalized Bolza problem of the classical calculus-of-variations. 

The powerful results associated with this classical theory thus are avail- 
able for attacking optimal control problems. 

The calculus-of-variations approach yields a nonlinear two-point 
boundary value problem for which closed-form solutions are usually not 
possible. Sophisticated numerical procedures have been developed because 
of the need to solve these problems. These procedures have become feas- 
ible in the last decade because of the development of large-scale digital 
computers . 

This chapter introduces the definitions and terminology used 
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throughout the dissertation, A brief history of the development of the 
numerical optimization methods is given also with interest centering on 
the second-order variational methods. The class of trajectory optimi- 
zation problems to be solved is stated, as well as the associated non- 
linear two-point boundary value problem obtained from using a calculus- 
of-variations approach . 

Chapter II discusses (1) the second-order variation approach 
used to solve the nonlinear two-point boundary value problem, (2) the 
general set of perturbation equations for second-order methods, and (3) 
techniques used to achieve an integration of these perturbation equations. 

Chapter III presents the Modified Sweep Method based on the 
generalized Riccati transformation. Chapter IV develops the linear feed- 
back control law for the Modified Sweep Method. The numerical results 
obtained using this new method are discussed in Chapter V, with conclusions 
and recommendations presented in Chapter VI. 

1.1 Definition of Terms Used 

The MSM (Modified Sweep Method) is obtained from the SSM (Suc^ 
cessive Sweep Method) by requiring that the control satisfy both local 
optimality and strengthened Legendre-Clebsch condition over the entire time 
interval of interest. This optimal control is then eliminated from the 
Hamiltonian for the problem and the restructured Hamiltonian used to obtain 
the nonlinear differential equations for both the state and Lagrange multi- 
plier variables. As in the case of the SSM, the MSM then uses the 
generalized Riccati transformation to solve the linearized two-point 
boundary value problem in terms of the state and Lagrange multiplier per- 
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turbations. This is done while simultaneously relaxing the terminal trans- 
versality requirements on the time-dependent Lagrange multipliers. It is 
desirable to compare the proposed MSM to other existing numerical opti- 
mization methods. For this reason, a study was made of several well-known 
methods which appear in the literature. This author felt that these methods 
were representative of the properties contained in the set of variational 
methods for the numerical solution of optimization problems. It is empha- 
sized that only a representative portion of the total number of existing 
numerical methods has been used for this study. In addition, the second- 
order methods intentionally have been selected more extensively than the 
first-order methods. The generalizations made, therefore, pertain only to 
those methods contained in Section 1.2 on the historical development of 
numerical optimization methods. 

This study of the selected group of existing methods revealed a 
set of properties which can be employed to describe the characteristic 
features for each method. These properties have been used to specify, ar- 
bitrarily, twelve classes of numerical optimization methods. Each method 
can then be identified as belonging to a particular class according to the 
following properties: 

(1) the ORDER (first or second) of the theory upon which the 

method is based. 

(2) the APPROACH (direct or indirect) used by the method to 

compute the required corrections. 

(3) the ITERATION PROCESS (interval, boundary or hybrid) used 


by the method. 
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The following definitions are used: 

Definition 1: Order of the Method 

A method is described as first order if it is based only upon 
the theory of the first variation for a real functional. If a method is 
based upon the theory of the second variation for a real functional, it is 
described as a second order method. 

Definition 2: Approach of the Method 

A method is said to take a direct approach if the required 
corrections (state, control or Lagrange multipliers) are computed such 
that the performance index for the augmented variational problem is 
itself directly affected in some manner to expedite convergence. If a 
method chooses to compute the required corrections based on the set of 
first-order necessary conditions required for optimality with respect 
to the control, then the method is said to take an indirect approach. 

Definition 3: Iteration Process for the Method 

A method is said to use an interval value iteration process if 
the end result of a particular iteration is the computation of correc- 
tions to the variables (state, control, or Lagrange multipliers) over the 
entire time interval of interest. Methods which compute corrections to 
these same variables at a boundary only are said to use a boundary value 
iteration process. If a method combines both an interval and boundary 
value iteration process, it is described as using a hybrid iteration pro- 
cess. 

Definition 4: Convergence for a Method 


(a) Control Function Iteration Method. Given an arbitrary 
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numerical tolerance e , a control function iteration method is said to 
have achieved convergence if 



where 

I Ih^I I = Max {Abs [H'*' (x,u,X ,t) ] } t < t < t, 
u(t) 

(b) Boundary Value Iteration Method . Given a numerical toler- 
ance e , a boundary value iteration method is said to have achieved con- 
vergence if 


1 1 E . 1 1 

+ 

| ImJ | 

+ £2? 

M -PM 


f 

f 


The symbols H, E^, M^ and £J^ are defined in Section 1.4. 

Using these definitions, Table I summarizes the twelve classes 
of numerical optimization methods extracted from the methods chosen as 
representative for the study. Reference numbers specify major studies in 
each class while the acronyms (SSM, etc.) identify the particular class 
for the three methods to be compared in detail. 
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TABLE I 

CLASSES OF NUMERICAL OPTIMIZATION METHODS 


CLASS 

ORDER 

APPROACH 

ITERATION PROCESS 

REFERENCES 

1 

1 

Direct 

Hybrid 


2 

1 

Direct 

Interval value 

4,12 

3 

1 

Direct 

Boundary value 


4 

1 

Indirect 

Hybrid 


5 

1 

Indirect 

Interval value 


6 

1 

Indirect 

Boundary value 


7 

2 

Direct 

Hybrid 


8 

2 

Direct 

Interval value 

19,26 

9 

2 

Direct 

Boundary value 

35 

10 

2 

Indirect 

Hyb rid 

3,13 

11 

2 

Indirect 

Interval value 

23,27 (SSM) 

12 

2 

Indirect 

Boundary value 

10,11,16 

(MSM,MPF) 


As shown in Table I, the class of second-order methods which 
take a direct approach in computing the desired corrections and implement 
a hybrid iteration process was not represented in the methods chosen for 
the study. Furthermore, Classes 1,3, 4, 5 and 6, all first-order methods, 
were also not represented. This can be attributed to the fact that second- 
order methods were of primary interest in this study. 
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The similarities and differences between the three numerical 
optimization methods (SSM, MSM and MPF) to be discussed are now obvious. 

All three are second-order methods which use an indirect approach in 
computing the required corrections for the state, control or Lagrange 
multiplier variables. The SSM falls into Class 11 because it uses an 
interval iteration process. Both the MSM and MPF fall into Class 12 
because each uses a boundary iteration process. 

A historical development of the methods chosen for the study 
is detailed now for reference purposes. 

1.2 Historical Information 

First-Order Methods . The first numerical procedure for solving 
control optimization problems which generated active interest was devel- 
oped independently by Kelley 12 and Bryson and Denham 4 . Their research 
extended the concept of steepest descent developed earlier by Courant 5 . 

The Class 2 method was based on the first-order variation of a scalar 
functional, with a control function assumed for the time interval of in- 
terest. Corrections to this control were then computed iteratively using 
an ordinary gradient technique. Applications showed that the method was 
easy to implement and tended to convergence with even gross initial control 
estimates. 

The method, however, possesses two undesirable features. First, 
the convergence rate decreased asymptotically during the terminal stages 
of iteration. Second, once convergence was achieved, the control obtained 
was only within a numerical tolerance of the Eulerian control. 

Due to the first undesirable feature, numerical procedures to 
increase the convergence rate flourished. These were all first-order 
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methods and are not discussed in this research. Both features led to the 
development of the second-order methods which sought to increase conver- 
gence as well as provide the true Eulerian control. 

Second-Order Methods . Jurovics and McIntyre 11 solved the two- 
point boundary value problem of trjaectory optimization by using the equa- 
tions which are adjoint to the linearized Euler-Lagrange equations. Their 
method was called the Adjoint Method (Method of Adjoint Functions), and 
extended the work of Goodman and Lance 7 to allow for variable terminal 
time. The Adjoint Method is a Class 12 method wherein an indirect approach 
is used to compute the desired corrections while iterating on the initial 
boundary values of the Lagrange multipliers. 

Breakwell, Speyer, and Bryson 3 developed a "second variation" 
method (Class 10) to solve control optimization problems. Kelley, Kopp, 
and Moyer 13 also developed a "second variation" method similar to the 
previous one. Jazwinski 10 developed a modified adjoint method equivalent 
to the second-variation method of Breakwell, Speyer, and Bryson 3 by ex- 
tending the method of Jurovics and McIntyre 11 . Jazwinski' s method had 
the specific advantage, however, of requiring considerably less storage. 
Furthermore, it required less computer time in that fewer integrations of 

an equivalent set of equations were necessary. 

McGill and Kenneth 20 developed the Generalized Newton-Raphson 

Operator Method for solving two-point boundary value problems. This method 
falls into Class 11, which uses the indirect approach linked with an 
interval value iteration process. A proof of quadratic convergence for 
the method had previously been given by these same authors 21 . This method's 
major drawback was the laborious manner in which corrections to the final 
time value were computed. 



An alternate approach, the Modified Generalized Newton-Raphson 
Method, using the same method was developed by Long 17 . To eliminate the 
awkward handling of free final time, a change of independent variable was 
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used for the free-final-time corrections. Another method based on the 
Newton-Raphson approach, but incorporating a better technique for compu- 
ting the free-final-time corrections, is the Modified Quasilinearization 
Method developed by Lewallen 16 . Sylvester and Meyer 36 have also used the 
Newton-Raphson approach, calling it quasilinearization. 

A method based on the theory of both the first and second vari- 
ations was devised by Merriam 26 . This is a Class 8 method in which a 
direct approach is taken for computing corrections to the control functions 
assumed throughout the time interval of interest. This particular method 
was instrumental in the development of the successive sweep method discus- 
sed in the next paragraph. 

McReynolds and Bryson 25 introduced the successive sweep method 
for solving optimal control problems. Although the method is a second- 
order method, the Eulerian control requirement is relaxed. The method is 
based on the generalized Riccati transformation and falls into Class 11 
(Table I) . A similar method called successive approximation was developed 
by Mitter 27 . He also showed the formal equivalence of this method to 
Newton's Method. 

Lewallen 16 also introduced the Method of Perturbation Functions 
(MPF) , based on previous work by Breakwell, et al. , 3 and Jazwinski 16 . The 
method falls into Class 12. Lastman 14 has shown the equivalence of all 
these methods to Newton's Method. 


Sutherland and Bohn 35 have recently developed a method which falls 
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into Class 9, which uses a direct approach to compute boundary corrections 
for the initial values of the multipliers. 

Mayne 19 developed a second-order method which falls into Class 8 
(Table I). However, a dynamic programming technique is used to attack the 
optimization problem instead of a calculus-of-variations technique. 

Jacobson 9 extended the Class 11 features into a new second-order 
algorithm through use of a differential dynamic programming technique. His 
method generalizes the successive sweep methods of McReynolds 23 and 
Mitter 27 . 

The development of the MSM completes the historical development 
for the numerical optimization methods chosen. 

As was mentioned previously, it is now desirable to compare the 
MSM to both the SSM and the MPF. Toward this end, a general class of con- 
trol optimization problems is first chosen. This class of problems is 
presented in the next section. 

1.3 Class of Control Optimization Problems to be Solved 

Posed as a special form of the Bolza problem from the calculus- 
of-variations (see Bliss 1 ), the general class of control optimization prob- 
lems to be solved is stated as follows: 

In the time interval t < t < t - , find an m-vector of control 
variables u(t) to minimize the real functional, 

t 

t 

o 


J(u) 


G(x f ,t f ) + 


Q(x,u,t) dt 


( 1 ) 
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subject to the n-vector of differential constraints. 


x = f(x,u,t) (2) 

while satisfying the p-vector of known initial conditions 

N(x ,t ) = 0 (3) 

o o 

and the q-vector of desired terminal conditions 

M(x f ,t f ) = 0 (4) 

The control and state variables in the following discussion are 
assumed to be defined on completely open regions and thus are not subject 
to inequality constraints. 

1.4 Associated Nonlinear Two-Point Boundary Value Problem 

Proceeding in the usual calculus -of-variat ions manner for solving 
the Bolza problem of control optimization stated in the previous section, an 
augmented functional denoted as I is first formed. This augmented func- 
tional has the property of being formally equivalent to the original func- 
tional; it incorporates the desired auxiliary conditions through the use 
of Lagrange multipliers. To form this augmented functional, the n-vector of 
Lagrange multipliers X(t) and the p and q vectors of constant Lagrange 
multipliers y and v adjoin the desired auxiliary conditions to the ori- 
ginal functional as follows: 


= J + 


T 

y N 


T 

v M 


X (f - x) dt 


I 


(5) 
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For convenience of notation, this functional is rewritten as 


I = P + ^ + j (H- Ax) dt 

t 


( 6 ) 


where 


P = P(x f ,v,t f ) = G(x f ,t f ) + v T M(x f ,t f ) 


V = V(x o ,n,t o ) = V T N(x o ,t o ) 


H = H(x,u,A,t) = Q(x,u,t) + A T f(x,u,t) 


and 


x = x(t) 


u = u(t) 


A = A(t) 


The scalar H is the variational Hamiltonian for this class of problems. 


Necessary Conditions . The set of first-order necessary condi- 
tions which must be satisfied by the extremal control for the augmented 
functional of the type above is obtained by requiring the first variation 
of this functional to vanish. These conditions are well documented in the 
literature (for example Bliss 1 , Hestenes 8 , Pontryagin 38 , et al^. , Tapley and 
Lewallen 38 ). In summary, these conditions are 


t 


x - H^(x,u,A,t) = 0 

A 


(7) 


t < t < t 
o - - 


A + H^(x,u,A,t) = 0 


H^(x,u,A,t) = 


( 8 ) 


0 


(9) 
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' - C o < 


E o ■ < r x + * T ) t ■ » 

o 


NCx t ) = 0 

o o 


fl o * <? t - H, t - 0 

o 


( 10 ) 

(11) 

( 12 ) 


t = t. 


E f 1 - »*>« = 0 


M(x f ,t f ) = 0 


fl f = (p t + H) t f * 0 


(13) 

(14) 

(15) 


Equations (7) through (9) constitute 2n+m Euler-Lagrange equations for 
this class of problems. Equations (11) and (14) are the p+q specified 
initial and final values of the problem state variables. The remaining 
equations form the 2n+2 set of classical transversality conditions from 
the calculus-of- variations . 

The control optimization problem thus is posed as a nonlinear 
two-point boundary value problem for the 2n+m variables x(t), u(t), 
and X(t) and the p+q+2 parameters y, v, t Q , and t^ in terms of 
the 2n differential equations (7) and (8) , the m algebraic equations 
(9), and the 2n+p+q+2 conditions in equations (10) through (15). 
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It is assumed that the initial time t and n values of the 

o 

initial state x ( t Q ) = X Q are specif ied. Equations (10) and (12) are 
then identically satisfied and therefore disregarded in subsequent dis- 
cussions. 

Sufficiency Conditions . To ensure that the control satisfying 
these first-order necessary conditions does indeed generate a weak mini- 
mizing solution, the second-order variation of the augmented functional 
I(u) must be positive everywhere in the interval of interest when it is 
evaluated along an extremal trajectory (Gelfand and Fomin 6 ). This re- 
quirement leads to the following additional set of second-order conditions 

1. Strengthened Legendre-Clebsch Condition 

The strengthened Legendre-Clebsch condition must be satis- 
fied everywhere in the interval of interest. Specifically, 
for any arbitrary change in control <5u(t) , 

6u T H 6u > 0 (16) 

uu 

is required. 

2. Jacobi (Mayer) Conjugate Point Condition 

The Jacobi (Mayer) conjugate point condition must be satis- 
fied everywhere in the interval of interest. This requires 
that no two points exist in the interval • t < t < t j which 
are conjugate to one another. 

The following restrictions on the definitions presented in Sec- 
tion I are subsequently assumed in this research: interval iteration 

corresponds to control function iteration and boundary iteration corres- 
ponds to Lagrange multiplier iteration. 


CHAPTER 2 


THE PERTURBATION EQUATIONS 

The second-order variational methods seek to solve the nonlinear 
two-point boundary value problem associated with trajectory optimization 
by solving an equivalent linearized problem in terms of perturbations in 
the problem variables. The six classes of second-order methods are dis- 
tinguished by approach and iteration process. Furthermore, each method 
for a given class is distinguished by the technique used to perform the 
numerical integration of the perturbation equations which are obtained 
from a first-order perturbation of the Euler-Lagrange equations, viz.. 
Equations (7) through (9). These perturbation equations can assume one of 
two forms with each form patterned by the iteration process selected for a 
given method. 

As in the case of the SSM (Successive Sweep Method) , if a con- 
trol function iteration process is used, a "PE" scheme is used to obtain 
the perturbation equations. The acronym "PE" designates the following 
procedure: perturb the Euler-Lagrange equations and then Eliminate the 

control perturbations. If a Lagrange multiplier boundary value iteration 
process is used, an "EP" scheme is used to obtain the perturbation equa- 
tions. The acronym "EP" similarly designates the following procedure: 
Eliminate the optimal control and then Perturb the resulting Euler-Lagrange 
equations. These two schemes for generating the perturbation equations 
for second-order methods are detailed below. 

2 • 1 The PE Scheme 

The PE Scheme is associated with control function iteration 


15 
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schemes such as the SSM. It proceeds by first perturbing the first- 

order necessary conditions for stationary control. Then the perturbations 

in the control are eliminated from the perturbed Euler-Lagrange equations 

for the state and Lagrange multipliers and the appropriate transversality 

conditions. This scheme tacitly assumes that the matrix H is nonsin- 

uu 

gular everywhere in the time interval of interest. The set of perturba- 
tion equations associated with this scheme is outlined in Appendix A and 
are summarized here as 


5x 


A 

N 

B 


f > 

6x 

+ 

r ) 

V 

<5X 


-C 

-A 1 


6X 


-w 

i J 

t 




w J 


L j 


( 17 ) 


where 


A = -H. H 1 H + H, 

AU UU UX AX 


B = 


-H H H , 
Au uu uA 


C = -H H 1 H + H 

XU UU UX XX 


V = 


-1 T 
H, H SH 

AU UU U 


W = 


-1 T 
H H 6H 
XU uu u 


and H = H(x,X,u,t) while the Hamiltonian partials , H^, etc., are 

evaluated on the known trajectory. The known trajectory is called a nomi- 
nal (or reference) trajectory. 

Note that for this scheme, the system of perturbation equations 
assumes the form of a set of inhomogeneous first-order linear differential 
equations with time-dependent coefficient matrices. 
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2.2 The EP Scheme 

The EP Scheme is used by the boundary value iteration methods. 

T 

It proceeds by first using the control optimality condition H u = 0 and 

T 

the strengthened Legendre-Clebsch condition 6u H 6u > 0 to eliminate 
the control from the Euler-Lagrange equations for the state and Lagrange 
multipliers, as well as the appropriate transversality conditions. The 
revised set of equations are perturbed then to obtain the following homo- 
geneous set of first-order linear differential equations. 


6x 


<5A 


-C 


-i T 


f \ 

6x 


6A 


(18) 


where 


A = H 


Ax 


B - H 


A A 


C 

H 


= H 


xx 


- H^x(t), A (t) , u(x,A,t),tj 


it 1 

and u(x,A,t) is the control obtained using the conditions that H u = 0 

and H >0. It is important to note that this scheme involves the as- 
uu 

sumption that the Hamiltonian H for the particular problem is structured 

such that H T = 0 and H >0 can be used to obtain the explicit rela- 
v u uu 


* T 

tion u(x,A,t). Implicitly, such a relation is assured if H 


u 


0 and 


H > 0 : however, such an explicit relation may be impossible to obtain 

uu 


for some problems. 
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2.3 Integrating the Set of Perturbation Equations 

It has been pointed out previously that second-order methods 
within a given class differ only in the technique that is used to inte- 
grate the set of perturbation equations. The two techniques presently 
available for accomplishing this integrating are detailed below. 

Explicit Integration . Methods using this technique choose to 
integrate directly the perturbation equations to obtain the perturbed 
values over the interval of interest. Some investigators have found (see 
for example the work of Merriam 26 ) that these methods suffer from numer- 
ical instabilities. Instability here is used in the sense that small 
errors in numerical precision will become exponentially very large over a 
long interval of numerical integration. The nature of these instabilities 
is associated with the numerical integration for coupled systems of linear 
differential equations which have split boundary conditions and admit both 
increasing and decreasing exponential solutions. 

Implicit Integration . Methods which presently use this technique 
are based upon a transformation process such as the generalized Riccati 
transformation. The technique consists of bypassing direct integration of 
the perturbation equations, and integrating a set of auxiliary variables. 
These in turn can be used to compute the perturbed values for the variables 
in the perturbation equations. Several advantages are claimed for this 
technique. First, the differential equations for the new auxiliary vari- 
ables have been reported to be more stable numerically than the original 
perturbation equations. Second, these auxiliary variables contain addi- 
tional intrinsic information about the optimal trajectory for the problem 
being solved. 
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A strong case concerning increased numerical stability in first- 
order linear two-point boundary value problems has been made by Rybicki 
and Usher 31 . However, work by Williamson 39 and this author has revealed 
that problems which have large differences in sign and magnitude for the 
eigenvalues of the coefficient matrix in the linear system of equations do 
not behave well numerically with the generalized Riccati transformation 
technique. This author’s opinion is that a valid generalized statement is 
yet to be made concerning the numerical stability properties of the Riccati 
transformation technique. 

That the auxiliary variables could contain additional intrinsic 
information certainly proves to be true. The earlier methods lacked in 
one respect: after convergence had been achieved, they required that post- 

convergence procedures be used to test the Legendre-Clebsch condition and/ 
or the Jacobi-Mayer conjugate point condition. Those methods which used 
the M EP" Scheme to generate the set of perturbation equations automatically 
took the Legendre-Clebsch condition into account when eliminating the con- 
trol from the set of Euler-Lagrange equations. However, the Jacobi-Mayer 
conjugate point condition still must be tested. This condition was often 
ignored and the converged solution was assumed to be a local optimum. 

However, using the generalized Riccati transformation technique 
on the perturbation equations provides the advantage of additional infor- 
mation for the current reference trajectory. Information is contained 
among these auxiliary Riccati variables for testing the Jacobi-Mayer con- 
jugate point and abnormality conditions from the calculus-of-variations 
(McReynolds 23 ) . It is well-known that the existence of a conjugate point 
precludes a trajectory from being optimal. The existence of such a con- 
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jugate point can be automatically detected continuously during the back- 
ward sweep process. This is accomplished by use of the fact that the 
matrix solution to the Riccati differential equation becomes unbounded at 
a conjugate point. 

On the other hand, the abnormality condition is equivalent to a 
certain matrix of these auxiliary Riccati transformation variables be- 
coming singular at the initial time. This information is important be- 
cause such a condition is tantamount to the inability in making correct- 
ions for values of the terminal constraints. This abnormality condition 
occurs for the Bolza problem if the boundary conditions at the final time 
are not linearly independent (McReynolds 23 ) . 

This leads to speculation concerning additional information 
about the reference trajectory which might be contained in the other Ric- 
cati transformation variables, either individually or in some combined 
form. To this author’s knowledge, little work has been done in attempting 
to extract such additional information. 

2.4 The Generalized Riccati Transformation Technique 

The generalized Riccati transformation is a transformation which 
changes the original two-point boundary value problem in terms of the 
coupled linear system of differential equations to an initial-value prob- 
lem having uncoupled variables and boundary conditions. This initial 
value problem is now stated in terms of n original problem variables , 
and in the general case, a total of jn(n+l) + q(q+l)j / 2 + f(nxq) + 

3(n+q) + 2j auxiliary Riccati variables where n is the number of state 
variables and q is the number of terminal constraint. Since the coupled 
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system of perturbation equations is integrated implicitly by integrating 
these auxiliary variables, it was expected that the differential equations 
for the auxiliary variables would be numerically more stable than the ori- 
ginal equations. Furthermore, Breakwell and Ho 2 have shown agreement with 
McReynolds 23 in that the conjugate point condition is related directly to 
the boundedness of an (nxn) matrix of Riccati variables which must satisfy 
a matrix version of the scalar Riccati equation over the time interval of 
interest. 

This transformation approach proceeds to solve the original non- 
linear two-point boundary value control optimization problem in the fol- 
lowing manner. A solution to the original nonlinear problem is assumed, 
and the corresponding terminal conditions are obtained. In general, these 
conditions are not satisfied to within the specified error tolerances. 
Desired changes, in these terminal conditions are specified, and the gen- 
eralized Riccati transformation is used then to generate a linearized field 
of solutions about this assumed solution. The transformation allows the 
specified changes in the set of terminal conditions to be mapped back to 
the initial time, when the particular member of the field that also satis- 
fies the initial conditions is selected. A new solution to the original 
nonlinear two-point boundary value problem is then computed using the lin- 
earized corrections, which are obtained through use of the auxiliary Ric- 
cati variables. As before, the new solution does not satisfy the desired 
terminal conditions exactly due to the linearity assumptions. However, the 
process can be applied iteratively until the desired terminal conditions 
are satisfied to within a suitable error tolerance. 

Historical Background . A Riccati transformation technique was 
first used by Gelfand and Fomin 6 in their successive sweep procedure of 
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solving two-point boundary value problems for linear inhomogeneous systems 
of second-order differential equations. The same transformation was gen- 
eralized and discussed for systems of first-order equations by Rybicki and 
Usher. 31 McReynolds 23 ,24 and Mitter 27 used the successive sweep method 
with the generalized Riccati transformation technique to solve the non- 
linear two-point boundary value problem of control optimization. Schley 

apd Lee 32 have developed a Newton-Raphson method which uses the Riccati 
transformation technique. Speyer and Byrson 34 have extended the concept 

of the Riccati variables for the case when some of these variables become 
unbounded. Narha and Berry 28 and Omicioli 29 have applied the successive 
sweep method of McReynolds to the shaping of optimal finite-thrust orbit 
transfer trajectories for which the control function is characterized by 
discontinuities. McGregor 22 has used the same method but has introduced 
modifications to handle problems with inequality constraints which contain 
the control explicitly. Most recently, Longmuir and Bohn 18 have shown 
how this technique can be used with any second-order method. 

Analytical Development . The generalized Riccati transformation 
for the linearized control optimization problem can be written in matrix 
form as 


6X(t)' 


«x(t )' 

dM f 

= R(t) 

dv 



dt f 

V 


+ P(t) 


( 19 ) 
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where dM^, dfi^, dv, and dt^ are constants for a particular iteration 
and 


f 

K 

D 

£ 1 


f \ 

T1 

E 

F 

g 

p & 

c 

T 

T 



♦ 

y 

z 

s 






w 


T 

where K, E, y respectively map given state perturbations Sx(t) into 
changes in the multipliers 6A(t), terminal state dissatisfactions dM^ , 
and terminal Hamiltonian transversality dissatisfaction d£2^ and 

K(t) is an nxn symmetric matrix 

E(t) is a qxn matrix 

y(t) is an n-vector 

T 

Also, D, F, z , respectively map changes in the multipliers dv into 
changes in the multipliers SA(t) , terminal state dissatisfactions dM^ 
and Hamiltonian terminal dissatisfaction dfi^ , and 

D(t) is an nxq matrix 

F(t) is a qxq symmetric matrix 

z(t) is a q-vector 

The scalars £, g and s respectively map changes in the final 
time dt^ into changes in the multipliers 6A(t), terminal state dissat- 
isfactions dMj and the dissatisfaction in the terminal Hamiltonian, dft^ 
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The quantities n , C and <t> respectively map changes in termi- 
nal local optimality dissatisfaction or terminal transversality dissatis- 
factions by the multipliers into changes <5A(t), dM f , and d£2^ . 

Differentiating the generalized Riccati transformation (Eq. 19) 
with respect to time gives 


* > 

6 A 


t > 

<5x 


* N 

fix 

0 

• 

* R 

dv 

4* R 

0 

0 


dtf 


0 

S. J 


( 20 ) 


Expanding the middle term on the right and transposing to the left yields 


I 

n 

0 

0 



f<$x I 


R 


dv 


+ P 



( 21 ) 


Using the perturbation equations, Equations (17), to eliminate 6 A and 
fix leads to the following expression 



-C 
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5A 

fix 


-w 


v 



6x 

k = R 

dv 


dtf 

v - 


+ P 


( 22 ) 
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Now 
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6 A 
6x 


-C ~A 
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r 'v 

Sx 


6 A 
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(23) 


Using the first row of the Riccati transformation. Equation (19), the 
following relation can be readily obtained. 


Sx 
S A 
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0 
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f \ 

Sx 


dv 

dt 


fJ 




Substituting Equations (23) and (24) back into Equation (22) gives 


(24) 
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5x 

dv 
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+ P 


(25) 
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Multiplying and collecting terms for arbitrary 5x, dv and dt^ , 
the following equations must hold 


R - -S WT 


where 


p = -Sr 


T 

e y 


0 0 


An + W 
Bn + v 


Performing the matrix multiplications yields the familiar set of equa- 
tions for the Riccati variables. 


K D 


(A T + KB) K + KA + C j (A T + KB) D | (A T + KB)& 
E(A + BK) | EBD | EBJl 

y T (A + BK) | yT B D I yTBfc 
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/• N 

n 


'(A T + KB)q + 

(Kv + w) 

d 

dt 

C 

= - 

E(Bn + 

v) 


<f» 


y T (Bn + 

v) 


From Equation (28) , the following rates of change for the Ric- 

cati variables are found to be equal: E^(t) = D(t), y(t) = j,(t) and 

z(t) = g(t). If, then, at the terminal boundary E^(t^) = D(t^), y(t^) = 

T 

£(t^) and z(t^) = g(t^), the following will be true: E (t) = D(t), 

y(t) = SL (t) and z(t) = g(t). This means not only that the matrix of 

* 

Riccati variables R(t) given in Equation (19) is symmetric but also that 
Equation (26) itself is also symmetric. In this case, 


• T 

R = -S WS 

where S and W are defined on page 26* 


(30) 


Terminal Boundary Conditions . The derivations of the general- 
ized set of terminal boundary conditions for the Riccati variables are 
presented in Appendix B. In summary, these boundary conditions are ob- 
tained from 


dH f 

dfl f 
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( 31 ) 
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The boundary conditions thus are 


R(t f ) 
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(33) 



29 


Z? + H a 

f Dt u uu uA f 
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H H 6H” 
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(x T - H H _1 H ,) dZ, 
u uu uA f 


Interval Value Process . . Methods using this process start by 
assuming a control function over the time interval of interest. In 
general then, H u (t) ^ 0. For such methods, the Lagrange multipliers 
can be made to satisfy the transversality conditions identically; i.e., 
Ej = 0. The Successive Sweep Method is an example of a control function 
iteration process. The terminal boundary conditions for the Riccati 
variables are then obtained from 
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where 


t = f-H H”^(H + H ,P )) 

2 ( u uu ux uA xx J 1 
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Note that the terminal matrix R(tj) of Riccati variables is 

symmetric only when the control function satisfies the local optimality 
T 

condition that H (t.) = 9. This same conclusion has been drawn by 
u f J 

McGregor 22 . However, this contradicts the results presented by McReynolds 
and Bryson 25 and Hitter 27 . This argument needs to be resolved. 

T 

Boundary Value Process . These methods assumed that H^(t) = 0 , 
where t < t < t^ . In this case, initial values are guessed for the 
Lagrange multipliers, and subsequent corrections are computed iteratively. 
In general then, + 0: likewise, dE f 4 0. The terminal boundary 

conditions in this case are obtained as 


( \ 
6 A,. 
r 


f(p ) 

XX f 

T 

M 

X f 

a f 


f N 

6 Xf 


f \ 

-dZ f 

dM f 


M 

X f 

0 



dv 

+ 

0 



T 

[°f 

•T 



d tf 

t j 


•T 

— x^dZ f 

< 4 


(35) 


Note then that there is a basic difference in philosophy between a control 

function and a boundary value iteration process. In the first case, the 

T 

optimality condition that = 0 is relaxed at each point in the inter- 

val to ensure satisfaction of the transversality condition = 0 by 

the Lagrange multipliers. In the second case, the optimality condition 
T 

H - 0 is satisfied at each point in the interval while the terminal 
u 

transversality - 0 is relaxed on the multipliers. 



CHAPTER 3 


THE MODIFIED SWEEP METHOD 

Merriam 26 and Mitter 27 have pointed out that boundary-condition 
iteration methods have certain programming advantages ; viz . , computer 
logic is relatively simple, and programming storage requirements are 
small. Furthermore, accurate trajectories are obtained in problems where 
these methods are successful. Experience has shown that such methods have 
one main disadvantage, viz . , the numerical instability mentioned previously. 
The nature of this instability has been discussed by several researchers, 
among them Rybicki and Usher. 31 Since the Riccati transformation tech- 
nique attempts to circumvent this problem by dealing with new uncoupled 
variables, this approach enhances the desirable features already known 
about boundary iteration schemes. 

3.1 Differential Equations 

For the modified sweep method, it is assumed that the Euler- 

Lagrange equations are satisfied over the time interval of interest. Fur- 

T 

thermore, the optimality condition that H^ * 0 is assumed to yield an 
explicit expression for the control in terms of the other variables. The 
Eegendre-Clebsch condition is then used to yield the extremal control 

& 

u = u(x,X,t) (36) 

This expression can now be used to eliminate the control from the original 
nonlinear Euler-Lagrange equations for x and X as well as from the 
appropriate transversality conditions. The set of first-order necessary 
conditions can now be rewritten as 


31 
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C o - C ^ h 


x - H^(x,X,t) = 0 

X+H*(x,X,t) = 0 


t = t 
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N(x ,t ) = 0 

o’ o 

T T 

K = (p - x 1 ) = o 

f x 

M(x f ,t f ) = 0 

= (P t + H) = 0 


(37) 

(38) 

(39) 

(40) 

(41) 

(42) 


Equations (37) and ( 38 ) are 2n equations for the 2n unknowns x(t) 
and X (t) . Equations (39) through ( 42 ) are 2n+q+l conditions for the 
2n unknowns x(t), X(t), and the q-KL unknown parameters v and t^« 
These equations constitute the familiar nonlinear two-point boundary value 
problem. A first-order perturbation of the nonlinear Euler-Lagrange equa- 
tions is now considered. This yields the following homogeneous linear 
system of equations (see Section 2.2). 
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As was mentioned on page 27 , the differential equation for the 

n+q+1 matrix of Riccati variables, R(t) , will be symmetric if the terminal 

T 

boundary values are such that R(t^) - R (t^) . In the next section it 
is shown that R(t^) will always be symmetric for the MSM. This 
reason, along with the fact that an EP scheme is used by the MSM to obtain 
the perturbation equations given in Equation (18) , give the following dif- 
ferential equations for the MSM Riccati variables: 
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where 
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(45) 


3.2 Boundary Conditions 

Boundary conditions for these equations are obtained by a first- 
order perturbation of equations (39) and (.40) 

6x(t Q ) = Sx o (46) 

S1 f * <P ** ) f 5x £ + M s £ dv + “f dt f - dE f (47) 

Equations (43) represent 2n equations for the 2n unknowns <5x(t) and 
6A(t). Furthermore, equations (46) and (47) give 2n split-boundary con- 
ditions for these variables in terms of the 2n known parameters 6x q 
and 6A^ plus the q+1 additional unknown parameters dv and dt^. The 
required additional q+1 conditions are obtained by also performing a 
first-order perturbation of equations (41) and (42) . It is shown in 
Appendix B that this procedure yields 
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In matrix form, these boundary conditions can be summarized as 
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(50) 


Note that this coefficient matrix is symmetric. Furthermore, if n values 
of the state are specified at t £ , E £ = dE £ - 0. For this case, the ter- 
minal boundary conditions reduce to the homogeneous form 
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To summarize. Equations (43), (46), and (50) constitute the linear first- 
order, two-point boundary value problem for the 2n functions <5x(t) , 
6X(t), and the q+1 parameters dv and dt^ in terms of the 2n+q+l 
specifiable parameters 6 x q , 6X^ , dM^ , and d$^. The computational 

procedure followed by the MSM will now be outlined. 

3.3 Computational Algorithm 

The modified sweep method can be implemented as follows: 


Step 1 - Assume n+q+1 values for X(t o ) » v and t^ . 

Step 2 - Integrate the nonlinear Euler-Lagrange Equations (37) and (38) 

forward from t to t c , viz., 
o f 


x = H^(x,X ,t) 


-H x (x,X,t) 


Step 3 - Test the error norm 
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If this criterion is satisfied, exit; otherwise, continue to 
the next step. 
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Step 4 


Set the terminal boundary conditions for the Riccati vari 
ables 
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* « 


•T 

-x 1 dZ- 

£ J 




tf 


(52) 


(53) 


Step 5 - Integrate the Riccati variables backward from the final to 

the initial time using the differential equations 

. T 

R = -SWS 

. X 

p = -Sr 

Step 6 - Compute the n+q+1 corrections 6A q , dv and dt £ using 

the generalized Riccati transformation (Equation 19), eval- 
uated at the initial time. These corrections are 
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dv - (F - gs 1 g T ) t * ( (dM^ - O -gs - d> ) 

o 

- (D T - gs 1 £ T ) 5x ) t 

o 

dt = s 1 ( (dfi, - 4>) -g T dv -& T 6x ) t 

1 o 0 

and 

SX * f K6x + Ddv + i dt, + n ) * 

° ^ f a 

where it has been assumed that 

dM f = -e M f 

•k 

dQj - -e 

* 

0 < e < 1 


Step 7 - Repeat from Step 2 using the new values 


x i+ 1 (t o ) 


X (t ) + 6X 

o o 


i+1 


= v + dv 


i+1 


t^ + dtj 


3.4 Computational Advantages 

The computational advantages of the MSM over the SSM are 
that it has to integrate n+q less variables and requires considerably 
less storage than the SSM. The exact comparisons are shown in Tables 
II and III. Table II shows the number of variables which must be inte- 
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grated by each method in forward and backward directions. Table III 
shows the storage requirements for each method. Only 3n + 2(q+l) quan- 
tities need to be stored in the MSM case. For the SSM, however, 
M(n(q+3)+m+n(n+l) /2) quantities have to be stored where M is the total 
number of points in the integration interval. A quick check for a 
typical reentry problem with M = 1,000, n = 6, q = 3, and m = 1 
shows the MSM requires storage of 26 quantities while the SSM must 
store 58,000 quantities. 

It is speculated that use of the SSM for large complex prob- 
lems such as the Apollo 3-D reentry will require fixed step-size integra- 
tion routines with a large enough step-size to remain within the computer 
storage limitations. Furthermore, the large step size may lead to unsat- 
isfactory numerical accuracy. 



39 


TABLE II 

MSM AND SSM VARIABLES TO BE NUMERICALLY INTEGRATED 


STANDARD SWEEP METHOD 

MODIFIED SWEEP METHOD 

Forward: x - n 

Forward: x - n 

Backward: X - n 

X - n 

„ n (n+1) 
K " 2 

Backward: K - ^ ^ 

D - nxq 

D - nxq 

Z - n 

£ - n 

n - n 

n - n 

F _ aia+il 
2 

_ q(q+l) 
2 

g - q 

g - q 

c - q 

c - q 

s - 1 

s - 1 

* - 1 

♦ " 1 

y - n 


z - q 


Totals : 

Totals : 

n(q+5) + (3q+2) 

n(q+4) + (2q+2) 

+ ~ (n(n+l) + q(q+l>) 

+ (n(n+l) + q(q+l)) 


Difference: (n+q) less variables 


Note: If all values of the term- 


inal state are constrained, i.e. , 


if q = n, then ri = 0, £ = 0, 


<(> » 0 and the difference increases 


to 2(n+q)+l less variables. 



TABLE III 

COMPUTER STORAGE REQUIREMENTS 


40 


STANDARD SWEEP METHOD 

At every point in the integration 
interval, the following values 
must be stored: 


x(t) 


n 

u(t) 


m 

x(t) 


n 

K(t) 


n (n+1) 
2 

D(t) 


nxq 

£(t) 


n 


Let M = total number of points in 
the integration, then 

M Jn(q+3) + m + 

quantities have to be stored. A 
typical reeptry problem has 

M = 0rder(l,000) , n = 6, q = 3, 
and m = 1. 

Thus, 58,000 quantities must be 
stored over the integration 
interval. 


MODIFIED SWEEP METHOD 

At the initial and final points only, 
the following values must be stored: 

x - n 
o 

X — n 
o 

v - q 



Only 3n+2(q+l) quantities need to 
be stored from iteration to iter- 
ation. 

Compare 26 quantities with 58,000 
for the Apollo reentry problem. 


CHAPTER 4 


THE MODIFIED SWEEP METHOD GUIDANCE SCHEME 

Initial conditions for dynamical processes are difficult to 
control in actual problems. Errors often occur which may be due to in- 
ternal mechanical causes such as premature cutoff by a thrusting rocket 
motor. Regardless of where these errors occur, they have the cumulative 
effect of causing a deviation from the intended optimal path. It is then 
desirable to use the known information about the path to recompute a new 
control program to accomplish the mission objectives. This is done by 
determining the control function corrections 6u as a function of the 
state perturbation; i.e., 6u = 5u(6x,t). This is a guidance problem in 
optimization theory. The guidance relations are now derived using the 
MSM. 

The MSM assumes that from the local optimality condition 

T T 

H u = 0 and the strengthened Legendre-Clebsch condition <5u H^Su > 0 , 

the minimizing control u(t) can be expressed as an explicit function of 

the state and Lagrange multiplier variables; i.e., 

u = U(x,A,t) (54) 

Perturbing this expression to first-order gives 

6u = U x 6x + U X <SA (55) 

The generalized Riccati transformation is 
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6X = K6x + 

D dv 

+ 

l 

dt 

f 


n 

(56) 

dM £ = D T 6 x + 
t 

F dv 

+ 

g 

dtf 

+ 

4 

(57) 

dJlf = &^<$x + 

T 

g dv 

+ 

s 

dt f 

+ 

♦ 

(58) 


Solving the last two equations simultaneously yields 


where 


and 


dv 


Trdx + Try + 'T z 
11 12 13 

dtf 

= 

TTfix + TTV + TTZ 
21 22 23 

y 

A 

(dM f - 0 

z 

A 

(dQ f - <f>) 

A 

A 

(F - gs 1 g T )~ 1 

7T 

11 

A 

A(gs"" 1 £. T - D T ) 

TT 

12 

A 

A 

TT 

13 

A 

-Ags 1 

TT 

21 

A 

-l f T , Tv 

-s (g IT + A ) 

IT 

22 

A 

-1 T 

-S g TT 

12 

TT 

23 

A 

-1/ T .. 

-s (g tt 13 -1) 


(59) 

(60) 


Using Equations (59) and (60) to eliminate dv and dt^ from Equation 
(56) gives 

SX = (K + Dir + Air )6x + (Dir + Jtir )y 

11 21 12 22 


+ (Dir + £ir )z 
13 23 


( 61 ) 
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Substituting this expression into Equation (55) then gives 

Su - fu + IL (K + Dir + ilir )}Sx 
1 x X n 21 J 

+ U. ((Dir + jin )y + (Dir + Jlir )z) (e>o\ 

* 12 22 13 23 v 1 

This equation represents the linear-feedback control law for continuously 
correcting the optimal control program for a given perturbation Sx(t) in 
the vehicle state. The new control program is then obtained by adding 
these corrections to the converged control program. Note that for a 
given perturbation 6x(t) , the coefficients are obtained by integrating 
the Riccati variables forward using the initial values corresponding to 
the converged solution. 

In guidance work it is desirable to know the control corrections 
in terms of a known time- dependent matrix and the initial vehicle state 
perturbations, i.e., 

Su(.t) = L(.t) <Sx(t Q ) ( 63 ) 

where L(t) is an explicit relationship between the Riccati variables 
from the converged optimal trajectory. Attempts to yield a relation such 
as Equation (63) were unsuccessful. Numerical studies using the MSM 
guidance scheme therefore were conducted with the assumption of a contin- 
uously correcting procedure. 

An immediate disadvantage is obvious in the guidance scheme 
represented by Equation (62). Since values for the converged Riccati 
variables are not stored by the MSM except at the initial point, these 
variables must again be integrated forward from this initial time re- 
gardless of when the perturbation Sx(t) occurs. A more detailed dis- 
cussion of this problem is contained in the next chapter. 



CHAPTER 5 


DISCUSSION OF NUMERICAL RESULTS 

The modified sweep method algorithm was programmed for the 
UNIVAC 1108 digital computer at the Manned Spacecraft Center in Houston, 
Texas. The integration schemes used follow. 

5.1 Numerical Integration Routines 

Fixed Step-Size Integration . Fixed step-size integrations were 
carried out using an Adams-Bashforth predictor-corrector procedure with a 
Runge-Kutta starter (Lastman and Fowler 15 ). The Adams predictor had a 
discretization error of o(h 5 ), and the Bashforth corrector discretization 
error was of o(h 6 ); h is the step-size. The Runge-Kutta starter had a 
discretization error of o(h 5 ). Partial double-precision arithmetic was 
used as follows: the values of the dependent variables were carried in 

full double precision, but the derivatives were evaluated and stored as 
single-precision numbers. This technique minimized the effect of round- 
off error. 

Variable Step-Size Integration . Variable step-size integrations 
were carried out using a predictor-corrector-starter procedure a$ mentioned 
above. However, the discretization error in all cases was of o(h 5 ). 

These integrations were carried out in full double precision (Schwausch 3 3 ) . 

5.2 A Brachistochrone Problem 

To compare Modified Sweep Method converged results with known 
analytical solutions for a problem of sufficient complexity, a class of 
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free-final-time Brachistochrone problems was posed as follows. 



Minimize the value of the final time t^ for a particle to fall 
along a frictionless path in a constant gravitational field from point 1 
to point 2 subject to the constraints 


where 


and 


X - 

V cos u 

II 

• >> 

V sin u 

x(t o ) - 

0.0 

rt 

O 

II 

1.0 


V = / 2g(y - a) 

3 = y o"2i 


The variational Hamiltonian is 


H = V(X cos u + X sin u). 
x y 


Two differ 
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ent cases were solved for the terminal constraint vector : 


and 


M* = x(tj-) - 5.0 - 0 

0 

0 . 


x(t^) - 5.0 
y(tj - 8.0 


The Modified Sweep Hamiltonian and its Partial Derivatives 

T 

Using the optimality condition H u = 0 and the strengthened 

T 

Legendre-Clebsch condition that Su H Su > 0 , the control vector can 

uu 

be eliminated to yield the following Euler-Lagrange equations 






where 


A £ A* + X2 

x y 


Furthermore, the second partial derivatives required by the perturbation 
equations are 


H™ 

xx 


t 

0 

0 


isi 

v 3 
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H_ 

Xx 




H 

XX 


<i> 

A3 




X X 
l y x 


x x , 

x y 


-X 2 

x 


Results 

The M* and M 2 cases were computed using the fixed step-size 
integrator mentioned on page 44 . For comparison purposes , these two 
cases were solved using the Method of Perturbation Functions program, 

MPF (see Lewallen 16 ). A step-size h = 0.01 second was used in all 
cases, with the initially assumed values of the unknown Lagrange multi- 
pliers and final time as follows : 


X° = -0.2368 sec/ft 

X° = -0.6095 sec/ft 

* = 0.5410 sec 

—8 

The convergence criterion e was specified as 0.1 x 10 . The correc- 

tion procedure used in all cases was 25%, 50%, 75%, and 100% from the 
fourth iteration onward. 

Robe o£_ Convergence . Figures 1 and 2 show plots of the 

terminal constraint norms versus time for the Brachistochrone problem. 
Both the MPF and MSM results are plotted. Figure 1 shows the error 



NORM OF THE TERMINAL CONSTRAINTS 
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ITERATION COUNT 


FIGURE 1 

M j CASE FOR BRACHISTOCHRONE PROBLEM 


NORM OF THE TERMINAL CONSTRAINTS 
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FIGURE 2 

M 2 CASE FOR BRACHISTOCHRONE PROBLEM 
f 
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norm for the case M^ = x(t f ) - 5.0. The two methods converge in seven 
iterations; however, the MSM consistently shows a smaller error than 
the MPF for each iteration. Note that the decrease in the error norm 
for the MSM is significantly less than the MPF for the last iteration. 

Figure 2 shows the error norms for the case M| = (x(t f ) - 5.0 
y(t^) - 8.0J ; the MSM error norm is not always less than the MPF 

error norm. However, the error difference is never large. The terminal 
stages of iteration reveal the same high rate of convergence as in the 
M* case. 

The MSM for this problem at its worst took 20% less compu- 
tational time. However, this figure is not considered significant be- 
cause two unrelated computer programs were used. 

Aoonraoy of Converged Results ■ > The MSM converged solutions 

gave six decimal place agreement for both the M^ case and the M^ case 

when they were compared to the known analytical solutions. For the M^ 

case, the initially guessed multipliers A° and A° were in error with 

the converged values 224% and 282%, respectively, the initial guess at 

the final time was 2% in error. For the M^ case, the initial guesses 

on the same A 0 and A 0 multipliers were 500% and 260%, respectively, 
x y 

In this case, the initial-guess error for the final time was 14%. These 
results are tabulated in Table IV. 

To summarize, the MSM has exhibited rapid terminal convergence 
and reasonable convergence envelopes for the free-final-time Brachisto- 
chrone problem. 
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TABLE IV 

MSM INITIAL-GUESS ERROR RESULTS 



M f 

Converged 

Values 

% Initial 
Error 


M 2 

f 

Converged 

Values 

% Initial 
Error 

x° 

-0.0689 

244% 


-0.0357 

500% 

X 






x° 

-0.1623 

282% 


-0.1726 

260% 

y 






tr. 

0.5271 

2% 


0.6277 

14% 

f 







The excellent results warranted further applications of the MSM 

to more complex problems whose analytical solutions were not known and 

which were of current interest to the space program. For these reasons, 

an Apollo three-dimensional reentry problem was chosen. 

* • 

5.3 Apollo Three-Dimensional Reentry Problem 

In the time interval t Q < t < t. , find the roll angle program 
g(t) which can be used to control an Apollo spacecraft so as to minimize 
the weighted sum of heating and acceleration effects 



/l 2 + D 2 + 

m 


_ 1 

X p ? V 3 
0 


I 


dt 
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where X q is a constant weighting parameter. Here, the first-term in the 
integrand serves to measure the acceleration effects due to aerodynamic 
forces while the second term measures the convective heating experienced 
by the spacecraft. This minimization is to be accomplished subject to the 
differential equations of motion given as follows 


/* • \ 

h 


( TT • ^1 

V sxn y 

• 

e 


V cos y cos A/ (R + h) cos A 

• 

A 


V cos y sin A/ (R + h) 

• 

V 


G sin Y - 5 

• 

Y 


(G cos y/V) + (V cos y/(R + h)) + (L cos B/V) 

• 

A 

i - 


(-V cos y cos A tan A/(R + h)} - ((L sin g/(V cos y)) 


The following initial conditions represent the reentry conditions 
for a space vehicle on an Apollo-type lunar return mission. 


o 

■p 

JrC 


'400,000 ft ' 


'75.757576 mi 

CD 1 
/-*N 

rt 

O 

V 


0.0° 


0.0 rad 

A(t o ) 


0.0° 


0.0 rad 

v(t o ) 


35,000 fps 


6.8181818 mi/sec 

Y(t o ) 


-6.5° 


-0.11344640 rad 

A ( t o) 


0.0° 

l J 


0.0 rad 

.. 4 


G = -p/(R + h) 2 

D = pSV 2 C D /2m 


where 
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L 

P 


pSV 2 C L /2m 


p o e 


* 

-3h 


and 


m = spacecraft mass (assumed constant) 


Optimal reentry trajectories were determined for two cases of 
terminal conditions 


M f 


h(t f ) - h f 


9(t f ) 


A(t f ) 


V(t f ) 


- e, 


- a. 


- V 


Y(0 ~ 


A(t f ) 


f 

A, 


and M| 


0(t f ) 

A(t f ) 

V(t f ) 

A(t f ) 


- 0 , 


- A, 


- V, 


- A, 


Values for the terminal state represent a typical set of conditions at 
drogue parachute deployment for the Apollo space vehicle 


/- \ 
h f 


'75,504 ft 

®f 


24.1° 

s f 


-0.6° 

*f 

= 

856 fps 

*f 


-44.3° 

4H 

l<J 

.. - -V 


-29.4° 

** 


Numerical values for the Apollo parameters were assumed as 
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S = 129.3 ft 2 = 0.46379993E-05 mi 2 

m = 204.0 slugs 

p = 0.27E-02 slug/ft 3 = 0.39743447E + 09 slug/mi 3 

*° -1 -1 
B = 0.42E-04 ft = 0.2217600E + 00 mi 

1/2 

X = 0. 24509804E-07 sec/ (slug-ftr' 

0 1/2 
= 0. 17809708E-05 sec/ (slug-mi) ' 

Vi = 0. 14076519E + 17 ft 3 /sec 2 - 0.95629856E + 05 mi 3 /sec 2 

R = 0. 20908800E + 08 ft = 0.39600000E + 04 mi 


Figure 3 shows the essential geometrical relationships between 
the state variables for the three-dimensional Apollo reentry problem. The 
variables chosen to specify the state of the point mass spacecraft were 
h = altitude, 0 = longitude, A = latitude, V = speed, y = angle of at- 
tack, and A = heading angle. The following assumptions have been made: 
the earth is a nonrotating homogeneous sphere with its center fixed in 
interial space. Furthermore, its gravitational potential is characterized 
by an inverse-square law and it possesses an exponential atmosphere. 


The Modified Sweep Method Reentry Hamiltonian . In Appendix C, 
the Apollo reentry optimization problem is restated. The mechanics of 
restructuring the problem Hamiltonian by use of local control optimality 
and the strengthened Legendre-Clebsch condition are shown. The Hamiltonian 
which is optimal with respect to the choice of roll angle 8 is given as 


follows : 
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H = (pSV 2 /c 2 + C 2 /2m) + X p 1/2 V 3 + X V sin y + 7 ^--. T 
L D 0 1 (R + h) 


( cos A 


(X - X sin A) + (X sin A + X ) | - 


l cos ^ 26 3 5 ) (R + h) 2 

X cos y 

t \ sin * + — — ) 


q 

- 7^ pSV fx C_V + — /x 2 + X 2 cos 2 yl 

2m [ 4 D cos y 6 5 ' J 


The Euler-Lagrange equations for this problem are then generated 
by taking first partial derivatives of H as follows: 



and 



where i,j = 1, 


6 . 


These results, along with the second partial derivatives H , H . and H, , 

XX XA A A 

which serve as coefficients for the matrix Riccati equation are also pre- 
sented in Appendix C. 

Results . Initial attempts to solve the Apollo three-dimensional 
reentry problem encountered some difficulties when the modified sweep 
method algorithm was used. Using the system of units ft/lb/sec, certain 
elements of the Riccati matrix K grew very large at the initial time. 
Because of these large values, the matrix F also became very large at 
the initial time. Consequently, when its inverse was used to compute 
initial-time corrections, they were so small that the initial values al- 
tered only in the seventh decimal place. As a result, the initial tra- 
jectory was essentially duplicated by subsequent sweeps. 
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An attempt to impose arbitrary bounds on these variables was 
tried for several bounding orders of magnitude ranging from 0.5 to 
1 x 10 6 . In all cases, every element in the Riccati matrix achieved the 
bounding value by the third sweep. An attempt was then made to generate 
more accuracy by cycling through the evaluation-correction procedure (EC) N 
of the fixed-step-size integrator several times. Values were tried for N 
ranging from 2 through 9. This effort to prevent the Riccati matrix 

from going onto the limiting boundary was unsuccessful. 

A scheme which used a scaled fractional part of the corrections 

<5X q was then attempted, and this did not eliminate the numerical diffi- 
culties with the Riccati matrix. The vector M f was then altered with 
respect to size and to choice of terminal state variables, neither of 
which was successful. The system of units then was altered to slug/mi/sec, 
for which the range of Lagrange multiplier magnitudes became smaller. Al- 
tering the unit system was tried after discussion with Williamson 39 , 
whose studies on the same problem with the MPF revealed a correlation 
between the numerical sensitivities of the Lagrange multiplier equations 
and the unit system chosen. The choice of slug/mi/sec achieved a more 
suitable scaling for the magnitudes of the multipliers; however, this did 
not succeed in eliminating the difficulties with the Riccati matrix. 

A variable step-size integrator routine was then introduced 
which revealed the numerical sensitivity of the Apollo three-dimensional 
reentry problem to the single-step error on the UNIVAC 1108. This sensi- 
tivity was measured by fixing the final time at t^ = 437.263 seconds; 
the initial values for the state and Lagrange multipliers were defined as 
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X q and A q as shown in Table V. The state and Lagrange multiplier dif- 
ferential equations were then integrated forward from t =0 seconds to 

o 

t f . Using the terminal values for the state and Lagrange multipliers, 

the process was reinitialized at t f and a backward integration carried 

out. The values obtained at t Q using this procedure were then compared 

for agreement with the defined values of X and A . Four cases were 

o o 

tested in which the single step error e was bounded: 


Case I 1.0 x 10 10 < e < 1.0 x lO -07 


Case II 1.0 x 10 12 < s < 1.0 x 10 09 


Case III 1.0 x 10 13 < e < 1.0 x 10~ 10 


Case IV 1.0 x 10~ 14 < e < 1.0 x 10 _11 


To numerically integrate forward from the initial to the final time and 
reinitialize and integrate backward to reproduce the initial values to 
eight decimal places, the single-step error e had to be bounded as 


1 x 10 ^ < e < 1 x 10 H 


When the error became less than 1 x 10 , the integration step size 

was doubled for the next step. If the error exceeded 1 x 10 the 

step size was halved; otherwise, the step size remained unchanged. These 
results are summarized in Table V where the bar under the digits denotes 
deviations from agreement with initially assumed values. 

Figures 4 through 9 give a particular set of numerical results 
for this problem. This set of results was essentially identical for both 
set of terminal conditions M* and m| . 
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TABLE V 

SINGLE-STEP ERROR TOLERANCE 
INITIAL VECTOR 


0.71895168 - 02 


0.75757576 + 02 

-0.72100154 + 01 


0 

0.81929784 + 01 

X = 
o 

0 

0.24206058 + 01 


0.68181818 + 01 

0.16453211 + 02 


-0.11344640 + 00 

0.35760665 + 01 


0 


NUMERICAL RESULTS 



Case I 

Case II 

Case III 

Case IV 

ERRMIN 

1.0D-10 

1.0D-12 

1.0D-13 

1.0D-14 

ERRMAX 

1.0D-07 

l.OD-09 

1.0D-10 

1.0D-11 

X h 

0.719.34619-02 

0.71895286-02 

0.71895265-02 

0.71895128-02 

x e 

-0.72100154+01 

-0.72100154+01 

-0.72100154+01 

-0.72100154+01 

X A 

0.81929769+01 

0.81929783+01 

0.81929784+01 

0.81929784+01 

X V 

0.24207.833+01 

0.24206029+01 

0.24206023+01 

0.24206059+01 

Xy 

0.16458225+02 

0.16453227+02 

0.16453225+02 

0.16453212+02 

*A 

0. 357502.46+01 

0.35760662+01 

0.35760665+01 

0.35760665+01 

h o 

0. 7576.6004+02 

0.75757226402 

0.75757596+02 

0.75757572+02 

0o 

-0.40744522-06 

-0.13346897-02 

-0.15685918-08 

-0.14599980-09 

Ao 

-0.65481932-06 

-0.13262769-08 

-0.18274255-08 

-0.10809805-09 

V o 

0.68182440+01 

0.68181886+01 

0.68181827+01 

0.68181819+01 

Y° 

-0.11346347+00 

-0.11344642+00 

-0.11344644+00 

-0.11344640+00 

A 

0.41275475-05 

0.20899670-02 

0.10682438-02 

0.73654477-09 

o 







6 


h, Ml 

V, Ml /SEC 
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FIGURE 5 

REENTRY LONGITUDE, LATITUDE, ANGLE OF ATTACK, 
AND HEADING ANGLE VERSUS TIME 
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SLUG-MI 

SEC 



FIGURE 8 

X Y AND 10X A VERSUS TIME 
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Figure 4 shows the altitude and scaled reentry speed histories 
that the Apollo spacecraft would follow during optimal reentry for both 
the M* and m| cases. Figure 5 shows the histories for the longi- 
tude, latitude, angle of attack and heading angle state variables during 
these optimal reentries. These two figures are of interest because they 
define the optimal state histories which the Apollo spacecraft should fly 
to achieve the specified terminal conditions while minimizing aerodynamic 
acceleration and convective heating. Figures 6 and 7 show the X^ 
and X^ multiplier histories, respectively. These are the Lagrange 
multipliers associated with the rates of change in altitude and latitude. 
They are included here to define the trends to be anticipated for the 
specified set of initial reentry conditions. 

The two Lagrange multipliers which are required to define the 
optimal reentry roll profile are shown in Figure 8. This particular fig- 
ure shows the X and X. histories where X and X. are associated 
y A y A 

with the rates of change in the reentry angle and heading angle, respec- 
tively. Figure 9 shows the reentry history for the payoff function. 

The aerodynamic acceleration and weighted convective heating experienced 
by the spacecraft have been plotted to reveal their individual character- 
istics. A study of this figure shows that two peaks occur in both space- 
craft acceleration and heating during the optimal reentries. The optimal 
reentry roll procedures seem to call for a trade-off philosophy between 
acceleration and heating experienced by the Apollo spacecraft. The high 
peak in convective heating is initially balanced with a smaller acceler- 
ation peak in the vicinity of 100 seconds. This situation is subsequently 
reversed in the vicinity of 400 seconds where the high acceleration peak 
is balanced with the smaller heating peak. 
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Figure 10 and 11 show the X^ histories for the and M 2 

cases, respectively. This is the Lagrange multiplier which is associated 
with the rate of change of speed for the Apollo spacecraft. These figures 
were included to reveal, for the specified initial reentry conditions, the 
sensitivities of this variable to a change in terminal conditions. It is 
interesting to note that both histories are similar with the major dif- 
ference arising beyond 400 seconds. 

Figures 12 and 13 show the optimal reentry roll programs for 

the m| and M^ 2 cases, respectively. These are the roll profiles that 
an astronaut would have to use during reentry from a lunar mission to 
minimize aerodynamic acceleration and convective heating while satisfying 
the desired terminal constraints on the vehicle state. 

In each case, optimal reentry calls for the spacecraft to com- 
mence the reentry maneuver with the lift vector pointed almost straight 
downward. The spacecraft is then quickly rolled such that the lift vector 
is pointed almost straight up after 90 seconds. A slower downward roll 
of the lift vector is then initiated so that this vector is at a value of 
169 degrees by 350 seconds. The lift vector is subsequently rolled up- 
ward to approximately 15 degrees by 410 seconds at which time terminal 
downward roll procedures differ depending on the specified values for the 
final vehicle state. Specifying two less conditions on the final state of 
the Apollo spacecraft calls for less terminal roll' of the lift vector. 

An attempt was made then to obtain a precise evaluation of the 
MSM computational characteristics. For comparison purposes, the 
three-dimensional Apollo reentry M 2 case was chosen. The MSM program 
was then altered to assume computational characteristics similar to 
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3 9 

Williamson s MPF program. Values for the initial Lagrange multipliers 
and final time were specified in varying degrees of accuracy ranging from 
four to eight decimal places. Both programs were run on the University 
of Texas CDC 6600 digital computer using single-precision arithmetic 
everywhere except in the variable step-size numerical integrators where 
partial double-precision was used. Each integrator had a single-step 
error tolerance of 1.0E-10< e< 1.0E-08. The correction procedure re- 
quired correcting 100% of the terminal error after each iteration. 

Results of this comparison study are summarized in Table VI. The 


TABLE VI 
MPF/MSM 

Convergence Characteristics 
for Case of the Apollo Reentry Problem 


Significant 
Digits for 

V t f 

Time to Converge 
CDC 6600 
(Seconds) 

% More Time 
Required by 
MSM 

Number of 
Corrections 
Required 


MPF 

MSM 


MPF 

MSM 

8 

33 

49 

49% 

1 

1 

6 

66 

106 

60% 

2 

2 

4 

165 

207 

25% 

5 

4 
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MPF program required less time to converge in each case. However, with 
decreasing accuracy in the initial guesses for A q and t^, the MSM 
revealed a tendency toward fewer required corrections and more competitive 
time-to-converge. Terminal error norms of both the MPF and MSM pro- 
grams for the case of four significant digits in A q and t^ are shown 
by Figure 14 • 

No direct computational comparisons with the SSM were avail- 
able. McGregor 22 used the SSM to converge the three-dimensional 
Apollo reentry problem for the case of terminally specified values for 
0, A, and V. This particular case was converged using a fixed step-size 
integration routine. However, as was discussed previously, the MSM re- 
quired a variable step-size numerical integration scheme to preserve the 
numerical integrity of the state and Lagrange multiplier equations. In 
addition, the MSM failed to converge this particular case of the three- 
dimensional Apollo reentry problem. This failure is currently under study 
by this author. Numerical comparison between the SSM and MPF for this 
particular case of the Apollo reentry problem can be found in the study by 
Tap ley, et_al. 37 

5.4 MSM Guidance Results 

The MSM guidance scheme was implemented for the three-dimen- 
sional Apollo reentry M 2 case. A 5% perturbation in altitude, speed, 
and angle of attack was initiated at t = 0 seconds to study initial 
reentry condition perturbation effects. A similar perturbation was initi- 
ated at t = 75 seconds to correspond to initial peaks in spacecraft 
heating and acceleration. In either case, it was assumed desirable to 
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correct so as to satisfy the originally specified terminal conditions. The 
guidance scheme then assumed that 

y(t) - z(t) =0 (t o < t < t f ) (64) 


Control corrections reduced to 

Su = (U x + U x (K + Du n + Air 21 ))«x ( 65 ) 

Numerical results revealed that the MSM guidance scheme failed 
to satisfy desired terminal conditions for specified vehicle state pertur- 
bations. Investigation revealed that numerical instabilities arising 
from attempts to forward- integrate the matrix Riccati equation were re- 
sponsible for compromising effective terminal guidance. Further study 
to suppress these instabilities is needed to achieve an effective MSM 
guidance scheme. 



CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 

A new second-order variational method (the Modified Sweep 
Method) was developed for solving the two-point boundary value problem 
of trajectory optimization. If differs from the original Successive 
Sweep Method in that the iteration process is now associated with modi- 
fying the initial values of the Lagrange multipliers instead of the 
control function over the time interval of interest. This approach re- 
quires considerably less computer storage and yields the Eulerian control. 
The new method was tested successfully on several classes of problems. 

The following conclusions were reached about the Modified Sweep Method: 

CONCLUSIONS 

1. The method has appeal for problems in which knowledge of 
the Eulerian control is critical. The MSM yields the Eulerian control 
over the entire time interval of interest upon convergence. 

2. Significantly less computer storage than the SSM was re- 
quired. Only 3n + 2(q+l) quantities were required by the MSM algorithm 
to compute the desired corrections from one iteration to the next. This 

is a desirable characteristic for larger-dimensional problems and small- 
storage computers. 

3. The numerical integration of a least n+q less variables 
than the SSM is required. This feature is desirable because less com- 
putation time is required. 
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4. Rapid terminal convergence characteristics of second-order 
numerical optimization methods is retained by the MSM. 

5. The conjugate-point test feature contained in the original 
SSM is also retained by the MSM. 

6. A numerical comparison of the MSM with the MPF for the 
class of free final-time Brachistochrone problems revealed that the MSM 
possesses acceptable convergence envelopes and competitive time-to- 
converge features . 

RECOMMENDATIONS 


1. The basic nature of the generalized Riccati transformation 
technique for solving the linear two'-point boundary value problem of con- 
trol optimization should be studied. It is possible that other equivalent 
combinations might possess a better structure for solving the two-point 
boundary value problem than the combination presently being used. 

2. Sensitivity of the MSM algorithm to classes of problems 
should be determined. This recommendation is made because of the method's 
failure to converge the three-dimensional Apollo reentry problem when ter- 
minal state values are specified for longitude, latitude and speed. 

3. The correction procedure used with the MSM should be opti- 
mized such, that the largest allowable correction is always attempted during 
a given iteration. This should be accomplished for the obvious reason of 
reducing computational costs by requiring fewer iterations. 

4. Properties of the other Riccati transformation variables 
and their relations to the reference trajectory should be studied. Cur- 
rently, only information about the Jacobi-Mayer conjugate point condition 
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and the abnormality condition is being extracted. This information is 
contained in only two matrices of the many used by the Riccati transfor- 
mation technique. 

5. The MSM algorithm should be extended to treat state as 
well as control inequality constraints. The need for this extension is 
obvious since most practical problems are subject to such constraints. 



APPENDICES 



APPENDIX A 


THE INHOMOGENEOUS SET OF PERTURBATION EQUATIONS 

Let x(t), u(t), and A(t) be functions associated with an 
extreme trajectory for the functional to be optimized. With the assump- 
tions made in Necessary Conditions , page 12 , the following Euler-Lagrange 
equations are necessarily satisfied: 

i t - - - 

x = H^(x,u,A,t) = f(x,u,t) (A.l) 

X = -H^(x,u,A,t) (A. 2) 

0 = H^(x,u,A,t) (A. 3) 

where ( ) indicates that the variables are to be evaluated on the ex- 
treme trajectory. 

Now assume a nearby trajectory characterized by the 2n+m func- 
tions x = x + 6x , u = u + 6u , and X = X + <5A. Substituting x, u, 
and X into Equations (A.l) through (A. 3) and expanding into a Taylor 
Series to first order about this nearby trajectory, the following equa- 
tions are obtained 


Sx = H Sx + H. 5u (A. 4) 

AX All 

SX = -H Sx - H Su - H ,5A (A. 5) 

xx xu xX ' 7 

<SH T = H Sx + H Su + H ,SA (A. 6) 

u ux uu uA 

where the partial derivatives of the Hamiltonian H are evaluated along 
the nearby trajectory. 
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Making the assumption that the (mxm) matrix H uu is non- 
singular, the control corrections can be obtained from Equation (A. 6) as 

Su = H* 1 |6H T - H 6x - H 6X| (A. 7) 

uu ^ u ux uX J 

Using this expression to eliminate <$u from Equations (A. 4) and (A. 5) 
then gives 


r 9 

6x 



A B 


f N 

6x 

+ 

f N 

V 

• 


T 





6X 


-C -A 

V J 


. 6X . 


-w 

v* / 


(A. 8) 


where 


A - -H, H -1 H + H, 

Xu UU UX AX 


B a 


-H, H _1 H , 
Xu uu uX 


C = -H H -1 H + fl 

XU uu ux XX 


V = 


-1 T 
H, H 6H 
Xu uu u 


w = 


-1 T 
H H 6H 
XU uu u 


Equation (A. 8) represents the inhomogeneous set of linear perturbation 
equations used by second-order variational methods. For computational 
purposes, the following is used: 

SH T (t) = -e H^t) , 0< e < 1 (A. 9) 

u uu u~ 

BOUNDARY CONDITIONS 

Boundary conditions for these equations are derived in Appendix 
B. They are summarized here as 
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5x(t ) = 6x =0 

o o 


and 


(SX(t-) = (P ), 6x £ + M dv + a.dt. + t, 

f xx f f x f f f 1 


where the vectors a and t are also defined in Appendix B. 


(A. 10) 


(A. 11) 



APPENDIX B 


FIRST-ORDER PERTURBATION OF TERMINAL CONDITIONS 


To allow for changes in the variable final time from iteration 

to iteration, the following linear approximation is used throughout this 

section. For an appropriate variable; e.g. , v^, assume that 

dv,. = 6 v £ + v- dtj. . 
f f f f 

The transversality conditions on the terminal values of the La- 
grange multipliers are expressed by the condition 



A first-order perturbation of this condition gives 


(B.l) 


dE f - <*»>£*<£ + <P xv ) f dv + < P *t>f dt r dX £ (B - 2) 

Replacing dx^ and dX^ using the linear approximation stated in the 
first paragraph above, and grouping terms yields 


DP 

= (P ) _6x c + (P ).dv + I — ^ - X I dt, - 8Xj. 
f XX f f xv f \ Dt L f f 


(B.3) 


Replacing X by use of X = -H and transposing the 8X^ term to the 
left gives 

,T 

T | UJ 

-L- M A " J- I — 

f XX f f 


DP 

= (P + M T dv + (-^ + H T |dt. - dZ^ (B.4) 


Dt 


The required terminal conditions on the state variables are 


given by 


M f = M(x f ,t f ) = 0 


CB.5) 
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Similarly, a first-order perturbation in the variables yields 

dM, = M 6 + MLdt. (B.6) 

f x f x f f f 

The transversality condition on the final value of the Hamilton- 
ian is given by 

n f = [P t + H(x,u,X,t)^ (B. 7) 

After using the linear assumptions stated in the first paragraph of this 
section, a first-order perturbation in the variables gives 


dfi- = I(P. + H ) ] -fix. + (P„ ) -dv + (H ) _6u. + (H,) r SX 
f tx x /J f f v tv f u f f I f f 


(P,, + H ) + (P^ + H )x + H u + H X 

tt t tx x u X 


if 


dt, (B. 8) 


Eliminating 6u^ using Equation (A. 7), fiX^ using Equation (B.4), and 
collecting terms gives 


dftj = 


P„ + H - H H _1 H + 
tx X u uu ux 


[H, - H H _1 H . ) P fix, 

\\ U UU uX J XX j£ I 


+ 

tv 


t ■ H H ^ .1 P ] dv 

\X U UU u\l XV Jf 


-£r(P + H) + (h - H H -1 H ) I + H T ) ]dt 

(Dt t lx uuuuX/lDt x j L f 


fn h 1 6h t 1 -[ 

1 U uu ujf V 


H - H H H I dZ. 
X u uu u X f 


(B.9) 


Substituting x = H , grouping terms again, and rewriting gives 
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= 


a T - H H -1 (H + H P ) 

U UU UX uA XX 


6x r 
f f 


((^-) - H H -1 H M T ] dv 
( \ Dt ) u uu uA x Jf 

(&*t 


DP 

v *T x *T T -1 

+ H) + x — ^ + x H - H H "H a 
Dt x u uu uA 


If dt f 


H H _1 6H T - x T dl £ + H H _1 H , dE £ 
U UU u f u uu uX f 




(B.10) 


where 


A 

a = 


& ♦ «*) 

\Dt x/t. 


Manipulating the coefficient for the dt^ term, it is possible to rewrite 


this coefficient as 


(£(£ * <) - 


E f nfr ~ H H ~ 1h i a |c 
f Dt u uu uA Jf 


In matrix form. Equations (B.4), (B.6), and (B.10) can be written as 


> 

r 

6A- 


f 


dM. 


f 


dfl f 

> 



(P ) f 

XX t 


M 

X, 


a f + x 2 


M 


•T 

M f +T 3 


a f 


3 f + t 4 


■v 

<$x f 


* ' 
T 1 


dv 

+ 

0 


dtf 

k J 


T 5 

V 


(B.ll) 


where 


A ( DI>X + h t ^ 

“E - Iw + \j t 




£ (i + A 
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t » -fH H _1 (H + H ,P )1 

2 [ u uu ux uX xx Jf 

t = -fH H _1 H ^M T 1 

3 { U UU UA xjf 

x = -fj£ ~ + H H _1 H a) 

4 ( f Dt u uu uA Jf 

and 

t 5 = (E H’' 1 6H T )-fx T - H H -1 H A d£ 

& { U UU ujf y u uu uxy f f 



APPENDIX C 


A MODIFIED SWEEP METHOD HAMILTONIAN WITH 
FIRST AND SECOND PARTIALS OF APOLLO 
3-D REENTRY PROBLEM 


Problem Statement 


Minimize the real functional 




J = 


/ L 2 + D 2 


m 


+ X p 1/2 V 3 
o 


dt 


subject to 


h 

e 


V sin y 
V 

(R + h) 


cos y cos A 
cos A 


A 

V 


(R + h") 11 sin A 

- sin y - — 

(R + h) 2 m 


Y = - 


cos Y 


(R + h) V 


V , L cos S 

cos Y + 


(R + h) 


m V 


^ V cos y cos A sin A _ L_ sin g 

(R + h) cos A mV cos y 


and satisfying the end conditions 


x(0) = C (a constant vector) 

o 


M(Xf,tf) =0 £ a (qxl) vector j 
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where 


p o e 


* 

-6h 


L - 2 pSV2C L 


D - |pSV2C D 


The variational Hamiltonian is 


H = p e" gh SV 2 / C 2 + C 2 + X p 1/,2 e“ e ^ 2 ^V 3 + X V sin y 
zmo LD oo l 


V cos y 
(R + h) 


fsTi <X 2~ x 6 sln + tt 3 aillA+ V 


(R + h) 2 


^ cos Y 

X 4 sln Y + * 5 -y- 1 - 


y=- P e' eh SV 
2m o 


X,C_V - C T fx, cos $ - X- — n - 
4 D L \ 5 6 cos y 


)) 


Partials of the Apollo 3-D Hamiltonian With 
Respect to the Roll Angle 


H 8 


3H 

98 


x 

2m 


P 




sin 6 


cos 8 
cos y 


= 0 


This implies that 


X 


5 


sin 8 + X_ 
6 


cos 8 
cos y 


= 0 
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or 


Thus 


We require 


Thus 


tan 8 = - 


A cos y 
5 


sin 8 = 


±A§' + \| cos 2 y 
6 ^ 


and 


cos 8 


-X 5 cos 8 


±/A^ + A*: cos 2 y 

6 3 


Sufficiency Condition 


H 


68 


— p e" 6h SVC T 
2m L 


X 5 cos 6 - 


X 6 sin 6 


cos y 


H 66 » 0 


for a local minimum 


, „ , sin 8 

X_ cos 6 “ X, 

5 6 cos y 


< 0 


or, using H_ = 0 

p 


cos y 


±A^ + X| cos^ y 

6 o 


< o 


7T 7T 

Requiring < Y < “ y~ yields 
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0 < cos y < 1 

Therefore, the "+" sign of the radical must be chosen for optimality in g. 
Finally , 

/a 2 + A 2 cos 2 y 
6 t> 

-A 5 cos y 
Af + A 2 cos 2 y 

6 o 

The Optimal Hamiltonian With Respect to g 
Eliminating sin g and cos g then yields 


sin g 


OPT 


cos g 


OPT 


H = p e~ gh SV 2 /C 2 + cl + A p^ 2 e &( h / 2 >v 3 
zm o 


L D o K o 


. , „ . . V cos y (cos A /, 

+ x i v Sln T + oTTh) [jSTi (x 2 ' \ stn 


4 ) 


+ (A. sin A + A_) 

3 o 


(R + h)' 


A u sin y + X. 


COS Y 

V . 


p e" eh SV 
2m K o 


j||f . Q 

A. C_V + — — — A 2 + A 2 cos 2 y 
4 D cos Y 6 5 
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First Partials of H With Respect to x and X 
The Euler-Lagrange equations for the MSM are obtained from the 
first partial derivatives of the Hamiltonian H. These equations are given 
below, where 


(m \ 


• / 

kJ 

» 

li 

•H 

* 


and i,j = 1, ... 6. 


H . f . 4 £ e- 6 h s - | X p 1 / 2 e - B < h «>v 3 

3h 2m L D 2 o K o 


V COS Y 
(R + h) 2 


cos A 


« (1 0 — sin A) + (^o sin A + ^c) 

cos A 2 6 '3 5' 


2 y 


(R + h)' 


X, sin y + 
4 ' 


X 5 cos y i 


V 


+ $ e" eh SV 

zm 


X C^V + A 2 + X 2 cos 2 y 

4 D cos Y 6 5 ' 


H - II - o 
H x 2 36 - 0 


H 


x. 


3H 

3A 


V COS Y 

(R + h) 


cos A 


cos 2 A 


(X„ sin A - X„) 


* 


- 15 = — e” 6h SV/c 2 + C 2 + 3X p 1/2 e" 6(h ^ 2) V 2 


H = — = — e SWCf- + + 3X p 

x. 9V m L D oo 

4 


+ X, sin y + — g . _X ._ 

(R + h) 


cos A 


(X„ - X sin A) 


cos A v 2 6 
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+ (X 3 sin A + A 5 ) 


Ac cos y 


(R + h) 2 V 2 


i * ( C. 

' 2A C^V + 


n p e S 
2m K o 


4-b* • — /x e + ^ cos2 y 


cos y 


~ 3H , „ V sin y 

H = — = A, V cos y “ / rt 
x 5 3Y 1 (R + h) 


cos A 


(A. - A c sin A) 


cos A 2 6 


+ (A 3 sin A + A 5 ) 


) 

J (R + 


(R + h) 2 


A^ cos y 


- Y sin Y ] 


1 „ _-0h.„ r 
- Tr~ p e SVC T 
2m o L 


A 2 sin y 


cos 2 y 'A 2 + cos 2 y 


H 


3H V cos y 
X g 3A (R + h) 


S | - ° j (A 0 - A c sin A) + A cos A 
cos A 2 6 3 


H, 



_3H 

13A-, 


= V sin y 


H, 



_3H 
A3 A 


V cos y cos A 
(R + h) cos A 


H, 


■ ©- 


V cos y sin A 

(R + h) 


H, 



_3H 
1 3 A , 


B sin y - p e‘ eh SV 2 C 

(R + h) 2 2m 0 D 
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V cos y y cos y 

l 9 V 

(R + h > (R + h) 2 V 

-Ip e 

-Kvc 

2m p o 

** -/x| + A 2 cos 2 y 


V cos y cos A sin A 

WJ 

(R + h) cos A 


* SVC 

1 -gh L 
2m P o 6 


COS Y A 2 + A 2 cos 2 Y 

D J 


H Matrix 

A A 


where 


H 


XX 


a = - 


o 

2m 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


C L SVe 


* 

-gh 


0 

0 

0 

0 

a 

b 


( x ! + °° s2 ->) 3/2 


b - 


C SVe 
Li 


* 

-gh 


( A 6 + X 5 COs2 ^) 3/2 


A 2 cos y 


X X cos y 
5 6 


c = - 


'0 

2m 


C L SVe 


* 

-6h 


U 2 + A 2 cos 2 y 


13/2 


A 2 cos y 
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H.. Matrix 
Xx 


H. = 0 

X 1 X 1 


H =0 

X 1 X 2 


H =0 

X l x 3 


H > x = sin y 
X 1 4 


H, = V cos y 
X l x 5 


H, = 0 

X 1 X 6 


H 


X 2 X 1 


V 

cog 

(R + h) 2 


cos A 
cos A 


H 


X 2 X 2 


0 


H 


X 2 x 3 


V 

(R + h) 


• cos y cos A sin A 
cos 2 A 


H 


X 2 x 4 


cos y cos A 
(R + h) cos A 
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V sin y cos A 
(R + h) cos A 


V cos Y sin A 
(R + h) cos A 


(R + h)‘ 


cos y sin A 


sin A 


(R + h) 


sin y sin A 


cos y cos A 


sin Y + B SV 2 C D e‘ 


(R + h) 


H. =0 
X 4 x 2 


0 
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H. 


X 4 x 4 


0 * 

-2- svc n r Bh 

m D 


H 


Vs 


(R + h)' 


cos y 


H = 0 

X * x 6 


H 


X 5 x l 


V cos y + 2y cos y 
(R + h) 2 (R + h) 3 v 


. p * 

* o 


+ ^ 2m C L SVe 


-$h 


X 5 cos Y 


Vx 2 + X 2 cos 2 
6 5 


H, - 0 

X 5 X 2 


H. = 0 
X 5 x 3 


H 


X 5 x 4 


cos y y_ 


cos Y _ ^0 -eh 
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C T Se 
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H 
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X 5 cos Y 

+ X 2 cos 2 y 
5 
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where 


T 2 


A 2 sin y 


[A 2 + A 2 cos 2 y 3/2 


6 5 


H 


Vl 


„ . . . * p C T SVe 

V cos y cos A sxn A , * o L 

+ g — 


* 

-6h 


(R + h) 2 cos A 


2m 


cos y Va 2 + A 2 cos 2 y 
6 5 


H, =0 
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H 


X 6 x 3 


H 


_V cos y cos A 


(R + h) 2 


X 6 x 4 
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p C Se 
O L 


* 

-eh 
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H 
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V6 
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H Matrix 
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+ e ~ e~* S 2X C_V + — A 2 + X 2 cos 2 y 

2m 4 D cos y 6 5 ' 
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6 5 


?S°2-L _ sin_i (x . t sln 4) + x cos A 
(R + h)2 cos A 2 6 3 


X 2 X 2 


X 2 X 3 


0 
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H =0 

X 2 X 6 


H 


X 3 X 1 


V cos y 
(R + h) 2 


cos A /, . , ^ \ 

(A sin A - A ) 

cos 2 A 2 6 


H =0 

X 3 X 2 


H 


X 3 X 3 


V cos y 
(R + h) 


cos A 


COS" 


;( 


A (1 + sin 2 A) - 2A sin A 
2 6 


H 


X 3 X 4 


cos Y 

(R + h) 


C - ° - S — A (A sin A - A ) 
cos 2 A 2 6 


H 


X 3 X 5 


V sin y 
(R + h) 


cos A * -x \ 

(A sxn A - A ) 

cos 2 A 2 6 


H 


X 3 X 6 


V cos y 
(R + h) 


S — - — (A sin A - A ) 
cos 2 A 2 6 


H 


X 4 X 1 


4p SV e 4h - \ tl p l/Ve 4 ' h/2 > 

o „ 2 oo 

in 


cos y 

(R + h ) 2 


COS f (*2 - *6 sin A > + CI3 Sin A + A 5 ) 


cos 


2v X 5 COS Y 


(R + h) 3 V 2 


* P n Se 

+ 3 — 

2m 


* 

-3h 


2X.CLV + — A 2 + A 2 cos 2 y 

4 D cos Y 6 5 ' 
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cos y cos A ,, . , s 

(TTw ^7 ( \ Sln 4 - V 

cos z A * 


p s AlS e 4h + 6 ~ l/2 Ve 4(h/2) 
0 0 0 


„ f X cos y i P r 

— - Se 6h X C 

L 1.\2 it 3 _ 4 u 


(R + h) 2 


A cos y - -SiS-JC =214 (X - x sln 4) 

1 (R + h) l“* 4 2 6 


\ X sin y 

+ (XsinA + Xj H § 


(R + h) 2 V 2 


£ - 6 \ 


Xr: sin y 


[cos 2 y /X 2 + X 2 cos 2 yJ 
v 6 5 J 


tt cos y sxn A ,, , . . , 

H = /— — r (X. - X. sin A) + X. cos A 

(R + h) cos A 2 6 3 


- siT1 1 — j (X - X R sin A) + (X sin A + X_) 

(R + h) 2 l COS A 2 6 3 5 


(R + h) 3 


X^ cos Y “ 


X 5 sin y 
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♦ K, — c l 

0 2m L 


Xg S in y 


[cos 2 y A| + x| cos 2 y J 


H =0 
X 5 X 2 


H 


X 5 X 3 


V sin y 
(R + h) 


cos A 
cos 2 A 


(X 2 sin A - X 6 ) 


~ , sin y 

H * X, cos y — /p . il\ 
x^x^ 1 ' (R + h) 


cos A « » , i \ 

(X 2 - X 6 sm A) 


cos 


+ (X, sin A + X ) - 


f X 5 sinY 


5 J (R + h) 2 


V 2 


D * 


X? sin y 


cos 2 y cos 2 Y 


~ , „ V cos Y 

H = -X,V sxn y — /p _rv\ 
x $ x^ 1 ' (R + h) 


cos A ,, , . 

7 (X, - X. sin A) 

cos A 2 6 


+ (X 3 sin A + X ) 


(R + h ) 2 


X k sin y + 


X- cos y X 


_ 1° S Ve _eh C t 

2m SVe C L T 4 


where 


cos y 




»! 


cos* 


3/2 


*6 


X 5 


cos 4 * 
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(1+2 tan 2 y) + X 2 sin 2 y 


V sin y 

(R + h) 


Sri (X 2 - X 6 sin A) + X 3 cos A 


V cos y 
(R + h) 2 


S * --' V (X, - X, sin A) + X, cos A 
cos A 2 6 3 


V cos y sin A ,, . , x 

(E + h) 77 (X 2 S “ 4 - X 6 > 

^ cos z A 


cos y sm A ,, , . .x . , 

rr ' 1 l 7 (X„ - X c sin A) + X 0 cos A 

(R + h) cos A 2 6 7 3 


V sin y sin A ,, , . , 

7 (X 0 - sm A) + X. cos A 

(R + h) cos A 2 6 ' 3 


V cos y cos A ,, , AX , j 

/T> . , < r (X - X. sin A) - X sin A 

(R + h) cos A 2 6 3 
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